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Abstract 

We use the extension of scaled particle theory (ESPT) presented in the accompanying paper 
[Stillinger et al. J. Chem. Phys. xxx, xxx (2007)] to calculate numerically pair correlation func- 
tion of the hard sphere fluid over the density range < pa 3 < 0.96. Comparison with computer 
simulation results reveals that the new theory is able to capture accurately the fluid's structure 
across the entire density range examined. The pressure predicted via the virial route is system- 
atically lower than simulation results, while that obtained using the compressibility route is lower 
than simulation predictions for pa 3 < 0.67 and higher than simulation predictions for pa 3 > 0.67. 
Numerical predictions are also presented for the surface tension and Tolman length of the hard 
sphere fluid. 



I. INTRODUCTION 



During roughly the last half century, the statistical mechanics of many-body systems has 
incorporated and has benefited from various forms of computer simulation. Indeed, with the 
remarkable advances in computer power that have occurred during that period, simulation 
has grown in the depth and breadth of its contributions to a point of near dominance over 
more analytical approaches. The hard sphere model provides a clear example, where its rich 
repertoire of both equilibrium and non-equilibrium phenomena are now better appreciated 
and understood as a result of a wide range of simulation activities ^ 2 ' 3 ' 4 ' 5 . 

Nevertheless, continuing to pursue the more analytical aspects of many-body statistical 
mechanics remains a valuable complementary mission. This and the immediately preceding 
companion papei^ (hereafter referred to as "I"), fall into this more analytical category. The 
objective has been to deepen understanding of the equilibrium properties of the classical hard 
sphere system. The underlying strategy involves an extension (ESPT) of the "scaled particle 
theory" (SPT) originally presented by Reiss, Frisch, and Lebowita^. The mathematical 
strategy of our approach has been described in "I". The present paper "II" provides details 
of a numerical investigation based on the mathematical formalism developed in "I" . One of 
our principal objectives is numerical prediction of the hard sphere pair correlation function 
g(r, p) over a substantial range in density p. 

The following Section II describes the coupled equations that must be self-consistently 
solved. At the center of this procedure is a non-linear integral equation for g(r, p). Successful 
solution of this equation at high density demands use of somewhat sophisticated techniques 
that are explained in Section II. Section III presents results not only for the basic quantity 
g(r,p), but also for the pressure and surface free-energy properties of the hard sphere sys- 
tem. This investigation has inevitably raised some important questions, particularly about 
what might be the most profitable directions for future extensions of the approach. These 
questions and our conclusions based on the results obtained thus far form the subjects for 
Section IV. 
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II. NUMERICAL PROCEDURE 



The goal of our numerical calculations is to obtain the hard sphere pair correlation func- 
tion g(r, p). In "I" , we introduced the overall general strategy behind this calculation. Here, 
we outline the solution procedure in greater detail. The computation can be divided into 
two major parts. The first involves the calculation of the isothermal reversible work W(X, p) 
required to create a single spherical cavity in a hard sphere fluid, where A is the cavity 
radius, and the second involves the calculation of the double or bi-spherical isothermal and 
reversible cavity formation work Wiir, A, p), where r is the separation between the two cav- 
ity centers. By calculating the work done in these two cases, the pair correlation function 
can be obtained using the following relation {(3 = l/fe^T): 

g(r, p) = exp[-/W 2 (r, 1, p) + 2(3W{1, p)) (1) 

To use the above relation, we note that the cavity work terms W(X,p) and W<i{r, X, p) 
must be evaluated at A = a, the hard sphere diameter (i.e., A = 1, using a as the natural 
scale for distance). Using the extended scaled particle theory (ESPT) of "I", the single 
cavity work is described by the following set of equations for different ranges of A: 

(3W(X,p)= - In ( 1 - 4vrpA 3 /3 ] , f < A < ~ j (2) 
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- In f 1 - 4vrpA 3 /3 + 2n 2 p 2 tf X dr 
+r 5 /12 



4AV/3 - A 2 r 3 



9(r,p)) A\<^<^ 1/2 ) (3) 

J(p)A 3 + K(p)X 2 + L(p)X + M(p) + N(p)/X , ( i < A ) ( 1 ) 
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Therefore, the calculation of W(X = 1, p) requires finding the coefficients J . . . N such that 
H^(A,p) is a smooth function of A. Among these, the coefficient J is related to the pressure 
P, K to the surface tension 7 and L to its curvature dependence. The fitting of these 
coefficients is performed at A = 3" 1 / 2 . In other words, we use the expressions given in the 
equations|HlandE]to calculate values of the coefficients that ensure a smooth transition across 
A = 3 -1 / 2 . This is accomplished by equating the work function value, along with its first 
and second A— derivatives, in the two regions (equations El and 0] ) providing constraints on 
three of the five coefficients. The other two necessary conditions are obtained by evaluating 
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the pressure P, and hence J(p) from the pair correlation function 

t/ \ 47T_„ A-ko Sti 2 p 2 ,„ . 

J(P) = -jP P = 3 + 9" ^' p) (5) 

and the following condition arising from the density of particle centers pg(r, p) at the surface 
of the A = 1 cavity 

(3W'(r,X = l) = 4irpg(l,p) (6) 

We note here that in this W(X,p) calculation, g(r, p) is required as input to evaluate the 
necessary conditions. We will further address this issue later on in this section. 

Having obtained W(X = l,p), we turn our attention to the bi-spherical cavity work 
W%{r, A, p). This calculation follows along similar lines as in the single cavity case. In "I", 
we developed expressions for the Wzir, A, p) for two regimes: one where the range of cavity 
sizes (A) limits the maximum possible number of clusters of spheres of unit diameter that can 
be found inside the bi-spherical cavity to 2, and the other where the maximum occupancy 
is greater than 2. We display them here in a different form. 

2 „2 rr 

r-2A 



m' 2 (r.\.p) = - In <j 1 — pV(r, A) + / r/n 



8 2 1 1 \ 27rV f r+2X 

—A 5 - -A 3 (r - rx) 2 + -A 2 (r - r^ 3 - — (r - r x ) 5 J n^(n) + — ^- y dr 

8 2 1 1 \ 1 

—A 5 - -A 3 (r - r 2 ) 2 - -A 2 (r - r 2 ) 3 + — (r - r 2 ) 5 Jr 2 ^(r 2 ) for n = 2 (7 

/3W 2 (r, A, p) = ^-V(r, A) + ^^4(r, A) + X(p, s)X + Y(p, s) + forn > 2 (8 

where V(r, A) and A(r, A) are the volume and surface area of the double cavity, and s is a 
transformation variable defined as s = r — 2A. Physically, s describes the distance between 
the nearest points on the surface of the two spherical cavities, when they are disconnected, 
and the corresponding distance between the interpenetrating surfaces, when they overlap. 
In the former case, s is positive and it is negative in the latter case. Using this variable 
allows us to circumvent singularities arising due to the disconnection of the two cavities of 
size A at r = 2A. 

The boundary between the two regimes, as discussed in "I", for r > (3 1 / 2 - l)/2 is at 
A = 1/2, and for < r < (3 1 / 2 - l)/2 is given by A = (l/2)(3 1 / 2 - r) + (l/8)(3 1 / 2 - r) _1 , 
where we have expressed distances in units of a. We assume that equation |H] holds for all r, 
provided A > 1/2. Furthermore, coefficients J(p) and K{p) are the same as those obtained 



in the single cavity W(X,p) calculation. Therefore, there are three coefficients that need 
to be determined, namely X(p,s), Y(p,s) and Z(p,s). These are determined by equating 
the value of the work function W 2) and its first and second derivatives with respect to A at 
A = 1/2. 

Along with the transformation variable s, another variable t = 2A — 1, a rescaled and 
shifted cavity size, is also used. In these variables, the fitting described above is done for 
each s value at t — 0. The pair correlation calculation in equation therefore requires that 
the coefficients X, Y and Z be determined for various values of the intercavity distance s at 
t — 1. The range of s in this calculation is — 1 < s < oo, which corresponds to 1 < r < oo. 

Once all eight coefficients J(p), K(p), L(p), M(p), N(p), X(p,s), Y(p,s), and Z(p,s) 
are computed, the calculation of g(r,p) follows from equation Q by setting A = 1 in the 
single cavity work W(X,p) calculation and t — 1 in the double cavity W2(s,t,p) calculation. 
However, as noted earlier in this section, both work calculations require g{r,p) as input(see 
equations El and Ej). To accomplish this pair correlation function calculation, a numerical 
solution technique is employed by using an estimate of g(r,p) as input, and then iteratively 
correcting it after calculating g(r,p) using the procedure discussed above. An arbitrarily 
small tolerance is chosen to determine whether the calculation has converged by comparing 
the the calculated and input g(r,p) functions. 

An appropriate numerical method to correct the input g(r,p) using the calculated values 
was determined after some trials on our part. The first attempted technique simply used a 
linear combination of the n th step input and calculated values in some arbitrary proportion 
to obtain the {n + l) th step input g(r,p). In other words, 

g%?(r, p)=xx g^ ut (r, p) + (1 - x) x gf^^r, p) (9) 

where x defines the proportion of calculated and input value in the corrected input g(r,p) 
for the next stage. This method is only effective at low densities (p < 0.4; here and in the 
remainder of the paper, dimensionless densities are obtained by scaling with a, the hard 
sphere diameter), regardless of the value assigned to x. The reason for this might be the 
high degree of non-linearity in the problem, which makes it difficult to achieve convergence 
using such a crude technique if the initial input g(r,p) is not an accurate enough estimate 
of the actual function. In this connection, we used, as two alternative input g(r, p)'s, a 
polynomial extrapolation of the converged pair correlation functions obtained for densities 



5 



lower than the one being computed, or the converged pair correlation function of a lower 
density case. These changes are able to improve the overall speed and convergence of the 
calculation. However, beyond a density p = 0.6, even these methods are unable to ensure 
convergence of the calculation. 

Satisfactory behavior in both convergence and speed was eventually obtained by using 
a generalized minimum residual method (GMRES) to solve the system of equations. The 
conventional GMRES method was developed for systems of linear equations 8 . It can be 
shown that for linear equations, the solution can be approximated by the n th iterate of the 
vector in a Krylov subspace (defined below) which minimizes the residual. Mathematically, 
for a system of equations of the form Ax = b, where x is the variable to be solved, the n th 
Krylov subspace is given by 

K n = span{b, Ab, A 2 b, A n_1 6} (10) 

where span{zi, Z2, ■ ■ • } refers to the set of linearly dependent vectors which can be con- 
structed from a vector set Zi, Z2, ■ ■ ■ The solution using the GMRES method is a vector 
in this subspace x n eK n such that the norm of residual Ax n — b is minimized as part of a 
least-squares problem. In order to solve this problem, we need to find the orthonormal basis 
vectors of K n . The component vectors in the above equation (b, Ab, . . ., A n ~ l b), being 
linearly dependent, are therefore not suitable to form this basis vector set. Instead, the 
Arnoldi process^ is used to construct the necessary orthonormal basis vectors i>i, . . ., v n of 
K n . This process also produces an upper Hessenberg matrix H n of dimensions n + 1 by n 
such that 

Av n = v n+1 H n (11) 

The orthogonality of the basis vectors leads to a simplification in the least-squares problem, 
or in other words the objective function is altered as 

min || Ax n — b ||= min || H n y n — e\ || (12) 

x„eK„ y n (R n 

where 

y n = x n /v n (13) 

and ei = (1,0,0,..., 0) is the first standard basis vector of W 1 . At each step in the solution 
procedure, GMRES solves for y n and then uses equation UHl to find x n . This process is 
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repeated with increasing value of n until the residual || H n y n — e± || is within an assigned 
tolerance value. 

This general approach, although originally developed for systems of linear equations, has 
been extended to handle nonlinear equations 10 . The main difference lies in the objective 
function. In solving a system of nonlinear equations F(x) = 0, the least squares problem 
at each step in the calculation is formulated as mm Snt K n \\ F(x n ) + F'(x n )s n || where s n is 
related to x n by the equation 

s n = s + x n (14) 

where Sq is a chosen parameter. The GMRES method can be shown to find a solution to a 
linear system of equations in at most n iterations 8 . For the nonlinear system of equations, 
however, no such guarantees can be made. As a result, a maximum number of iterations 
(n max ) is set. If the residual is not within the set tolerance value after n max iterations, 
the calculation is restarted such that sq = s nmax . A maximum number of restarts is also 
used as input to the method so that the calculation can be terminated in the event that 
convergence is not achieved. The maximum distance r max was set at r = 9 for these g{r,p) 
calculations. The computation interval 1 < r < r max was divided into 2000 discrete intervals, 
and Simpson's rule was used in all integrations. 

The major thrust of this paper has been to calculate the pair correlation function g(r,p) 
for a hard sphere fluid based on the extended scaled particle theory (ESPT) of "I" . Using the 
methods discussed in this and the previous section, g(r,p) has been computed for different 
densities. As an initial estimate for each g(r, p) calculation, the converged g(r, p) for a lower 
density case was used. A tolerance limit of 10 -7 was set for each calculation to determine 
whether the residual was low enough. 

III. RESULTS 

The solution scheme described in the previous section was used to obtain the equilibrium 
pair correlation function g(r, p) for the hard sphere fluid. We compare the g(r,p) obtained 
for different densities in Figs. ^ El and H3 The same figures also compare the calculated 
values of g(r, p) with those obtained from molecular dynamics simulations 1 - 1 ' . In general, the 
ESPT prediction for g(r, p) is in very good agreement with molecular dynamics calculations, 
indicating the success of the theory in capturing fluid structure and its density dependence. 
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However, for p > 0.4, the theory underpredicts g(l,p), as will be discussed in connection 
with equation of state predictions. 

The computation time required for each successive calculation to converge increases with 
density. In this work, we have been able to continue the pair correlation calculations to 
densities that are in the metastable region of the accepted hard sphere phase diagram. The 
generally accepted freezing density for the hard sphere fluid is Pf re ezin g = 0.943 4 . In Fig. 
gcj the pair correlation function based on the ESPT for a density p = 0.96 is shown. The 
second peak appears to be almost triangular in shape instead of the smooth peak observed 
at lower densities. Note that the development of this near-triangular peak shape is also 
evident in the molecular dynamics simulation results. Calculations of the pair correlation 
function for p > 0.96 using the above discussed numerical scheme do not converge. Further 
refinement of the numerical approach is under investigation. 

The contact value g(u,p) can be used to calculate the pressure of the fluid P using the 
virial theorem, 



In Fig 01 the pressure is plotted against density for the ESPT along with predictions from 
other equations of state (EOS's) for the hard sphere fluid. The pressure prediction from the 
ESPT of "I" corresponds closely to molecular simulation results^ at low densities. However 
at higher densities, there is significant deviation between the two. In fact, we find that the 
ESPT pressure is very close to the corresponding quantity computed via the Percus-Yevick 
closure^ (pressure equation). This, in itself, is an interesting result in light of the fact — 
that the original SPT of Reiss et al. 7 produces the same EOS as the Percus-Yevick closure 
(compressibility equation)^. Both versions of SPT are unable to accurately predict the 
pressure at higher densities. The original SPT over-estimates the pressure, while the ESPT 
underestimates it. 

An alternative calculation can be conducted using the isothermal compressibility equation 



The derivative (dP/dp)x can then be integrated numerically from the ideal gas limit, \(3P = 
p; g{r,p) — 1 for r > 1, as p — > 0]. We use density steps of 5p = .02 in this calculation, 
ensuring that computed dimensionless pressure values, Pa 3 / (/c^T), are accurate within 10~ 2 . 



pP = p+(2/3)7ip 2 a 3 g(a,p) 



(15) 




(16) 
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A cut-off distance in r is selected such that contributions to the integral in equation 1161 are 
considered only when \g(r, p) — 1| at any peak of g(r, p) is greater than 2 x 10~ 3 . The 
pressure obtained by this route is shown in Fig. |5J In the figure, we compare the EOS 
predictions from equations ITol and ITol with the corresponding original SPT prediction 7 . At 
high densities, the ESPT calculation via the compressibility equation closely matches the 
original SPT, or equivalentlyi^, the Percus-Yevick compressibility calculation. However, at 
low density, there is significant difference between the two, with the ESPT compressibility 
calculation yielding lower pressures. The same figure with the accompanying inset also 
displays EOS values as obtained from molecular simulation results. The EOS prediction 
from ESPT through the compressibility route underpredicts the pressure at lower densities 
(p < 0.67) and overpredicts this quantity at higher densities (p > 0.67), when compared to 
simulation results. 

The surface tension of the hard sphere fluid can be estimated from the parameter K in 
the asymptotic form of the cavity formation work /3W(\,p) in equation H] 

where 7 is the surface tension. A comparison of 7 from the ESPT with that of the original 
SPT is shown in Fig. El The original SPT has been shown to predict, with reasonable ac- 
curacy, the surface tension of the hard sphere fluid, when compared to simulation results^. 
Both theories predict negative surface tension for the the hard sphere fluid with 7 increasing 
in magnitude with increasing density. Similar values of 7 are predicted by both theories at 
low and high densities. However, for intermediate densities, the ESPT predicts a smaller 
magnitude for the surface tension than the original SPT. Comparison with simulation re- 
su ltgl£*iLi2i i s shown in Fig. Here, the pressure prediction for the different theories and 
simulations is plotted on the abscissa. We observe that the ESPT is able to provide a more 
accurate prediction of the surface tension at low density and pressure. However as the den- 
sity increases, the original SPT is able to more accurately capture the dependence of surface 
tension on pressure. Equation ITol provides an alternate route for the calculation of the EOS 
using the ESPT. As shown in Fig. |H1 use of the compressibility route enables the ESPT to 
provide an improved representation of simulation results for the variation of surface tension 
with pressure than the original SPT. 

The Tolman length 5 is a measure of the dependence of surface tension on curvature. It 
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can be estimated from parameter L in the asymptotic equation 0] from the following relation 



5 = -^ (18) 

87T7 v ; 

For a hard sphere fluid 5 > 0. The predicted variation of 5 with density is shown in Fig. El 
The Tolman length increases monotonically upon compression. Comparing these values with 
those from the original SPT, we find that ESPT predicts a more pronounced dependence of 
the surface tension on curvature. 



IV. CONCLUSIONS 

A new scaled particle theory of hard sphere pairs& has been used to calculated the pair 
correlation function g(r, p) of the hard sphere fluid over a wide range of densities, covering the 
generally accepted stable and slightly metastable fluid states. Comparison with results from 
molecular dynamics simulations^ shows that the extended scaled particle theory (ESPT) 
is able to accurately describe the fluid's structure, but it underestimates the g(r, p) contact 
value when p > 0.4. The theory is able to accurately reproduce the shape of the peak 
corresponding to second neighbor shell at r = 2 which becomes sharper as the density of 
the hard sphere fluid increases, in agreement with simulation results. 

The ESPT has been used to compute the equation of state of the hard sphere fluid by 
using the virial theorem. The calculated EOS values correspond closely at lower densities 
with those obtained from molecular simulations^. However, there are significant differences 
between the two at higher densities. We find that ESPT consistently underpredicts EOS 
values from molecular simulation. This can be contrasted against the overprediction of the 
EOS by the original SPT. Comparing the predictions of ESPT with those from other theories, 
we find that the calculated EOS values closely match those computed using the Percus- 
Yevick approximation (pressure equation). We have also calculated the hard sphere equation 
of state using ESPT with the compressibility equation. The resulting prediction is close 
to the equation of state obtained using the Percus-Yevick closure and the compressibility 
equation, especially at high densities. 

The surface tension and Tolman length have also been calculated. As expected, ESPT 
predicts a negative surface tension for the hard sphere fluid. ESPT is able to more closely 
capture the pressure dependence of the surface tension at moderate pressures than the 
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original SPT 7 . However, at higher pressures, the original SPT provides better prediction. 
As discussed above, employing the compressibility route improves the EOS prediction of 
ESPT, leading to a closer correspondence with 7 vs. P surface tension predictions from 
molecular simulation. ESPT predicts a larger Tolman length than the original SPT, or 
in other words, the new theory predicts a more pronounced effect of curvature on surface 
tension than the original theory. 

Having studied a small portion of the metastable branch of the equation of state, future 
research directions might include an examination of the high density range of the phase 
diagram in an attempt to calculate the pair correlation function of both the hard sphere 
solid and the metastable fluid. Using solid and fluid forms of g(r, p) as initial inputs to the 
calculation is a possible strategy that will be investigated as a means of testing the present 
theory's ability to predict a freezing transition. Further numerical work is needed to ensure 
convergence along the high-density metastable branch. 

Another avenue of investigation is to extend the current hard sphere theory to study 
solvation effects. The new theory would include the effect of orientation of and attraction 
between interacting particles. There have been other attempts at developing such theories 
in the pas t 19 ' 20 . These have been successful in highlighting important phenomena associated 
with hydrophobic hydration or solvation. 

Finally, we point out that additional consistency conditions can be imposed in the itera- 
tion scheme described in Section III. The extent to which such additional constraints improve 
the accuracy of ESPT predictions of thermodynamic quantities such as the pressure and the 
surface tension deserves further investigation. 
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FIG. 1: The pair correlation function obtained using the extended scaled particle theory (ESPT). 
Here, three densities are shown: (a) g(r, p = 0.2), (b) g(r, p = 0.3) (c) g{r,p = 0.4). The lines are 
ESPT calculations. The open circles are molecular dynamics simulations 11 . Distances have been 
scaled with the hard sphere diameter a. 
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FIG. 2: The pair correlation function obtained using the extended scaled particle theory (ESPT). 
Here, three densities are shown: (a) g(r, p = 0.5), (b) g(r, p = 0.6) (c) g{r,p = 0.7). The lines are 
ESPT calculations. The open circles are molecular dynamics simulations 11 . Distances have been 
scaled with the hard sphere diameter a. 
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FIG. 3: The pair correlation function obtained using the extended scaled particle theory (ESPT). 
Here, three densities are shown: (a) g(r,p = 0.8), (b) g(r,p = 0.9) (c) g(r, p = 0.96). The lines are 
ESPT calculations. The open circles are molecular dynamics simulations 11 . Distances have been 
scaled with the hard sphere diameter a. 
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FIG. 4: Equation of state of the hard sphere fluid. We compare predictions from ESPT using 
equation fTol (o). the original SPT of Reiss et al£ (0), the Percus-Yevick^ pressure (solid) and 
compressibility (dashed) calculations, and molecular dynamics simulation results^ (dot-dashed). 




FIG. 5: Equation of state of the hard sphere fluid We compare predictions from ESPT using the 
virial (o) and the compressibility (solid) routes, the original SPT of Reiss et alA (*) and molecular 
dynamics simulation results 12 (dashed). Inset shows details in the density range 0.4 < per 3 < 0.6. 
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FIG. 6: Surface tension of the hard sphere fluid. The solid line represents the prediction from 
ESPT; the dotted line, the original SPT of Reiss et ali£. 
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FIG. 7: Surface tension of the hard sphere fluid. The solid line represents the prediction from 
ESPT using the virial route; the dotted line, the original SPT of Reiss et alZ. Open circles (o) are 
simulation result o 16 i 17 i 18 . 
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FIG. 8: Surface tension of the hard sphere fluid. The solid line represents the prediction from 
ESPT using the compressibility route; the dotted line, the original SPT of Reiss et alii. Open 
circles (o) are simulation result o 16 i 17 i 18 . 




FIG. 9: Tolman length of the hard sphere fluid as predicted by ESPT (solid) and the orginal SPT 
of Reiss et al. 7 (dashed). 
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